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Abstract 


A mathematical framework is developed to solve the commercial aircraft trajectory optimization 
problem in a vertical plane for a general case. The generality refers to an arbitrary objective 
function, a spatial wind field, a dual-element control vector, and a set of constraints. The original 
problem is singular on one of the control elements and regular on the other. Due to the small 
aerodynamic path angle, which is typical of commercial aircraft flights, and a large time scale 
separation between the flight dynamics, it is conventional to take the aerodynamic path angle as 
a singular control. This leads to a model that is singular on both control elements. Followed by a 
theorem, we present a compact closed-form solution for the multi-control singular model in the form 
of Jacobian compositions through the Pontryagin Maximum Principle. The optimal structure is 
dependent on the boundaries as well as the relaxed parameters involved. An algorithm is proposed 
to capture various combinations of bang and singular extremals for the control elements, leading to 
some new optimal structures, as will be shown through the case study. Results due to the singular 
model are compared with the solutions due to the original model with the overall conclusion of 
extreme accuracy. Finally, the paper terminates with a discussion on the merits of the singular 
solutions in terms of computational time and accuracy. 


Keywords: Indirect Optimization; Optimal Control; Aircraft Trajectory Optimization; Closed 
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1 Introduction 


he aerospace industries are inquisitively in demand of modern optimization tools to minimize 
the total flight cost. The cost is usually a combination of the arrival time, fuel burn together 
with some environmental hazards such as climate and noise. 

Aircraft trajectory optimization as an Optimal Control Problem (OCP) has been the target 
of many research centers. The solution for this class of OCPs has been sought through the 
well-posed optimization methods; i.e., mathematical parameter optimization, heuristic parameter 
optimization, Pontyagin Maximum Principle, PMP, and dynamic programming (in approximate 
or classic form). Let us refer the mathematical parameter optimization as the direct mathematical 
(gradient-based) method or in short, the direct method and PMP as the indirect method, being 
two counterpart mathematical methods in optimization theory. 
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Due to the complexity and sensitivity of the resulting multiplier system from PMP, the indirect 
method shares the least contribution among the rest in aircraft trajectory optimization (see [1]). 
This complexity, by itself, has drawn researchers to the direct methods as well as the heuristic 
methods with limited level of accuracy and warrant. 

The direct method usually prevails in the matter of labor with a competing level of accuracy 
compared to its counterpart indirect method. The multipliers in a direct method imitate the 
presence of co-state dynamics through Karush-Kuhn-Tucker (KKT) conditions, which explains 
the accuracy of direct methods. Therefore, mathematically speaking, the direct method is a 
discrete version of PMP in such a way that the solutions from the direct method correspond 
to counterpart indirect solutions and the correspondence is regarded as the Covector Mapping 
Principle (CMP), a lemma, initially postulated by Ross ( [2], [3]). For this reason, the direct 
solutions are sometimes regarded as the numerical solutions, whereas the indirect solutions are 
the analytic ones, the terminologies used throughout this paper as well. 

The goal of the present research is to develop a mathematical framework to solve a general 
version of the commercial aircraft trajectory optimization problem in a vertical plane. Relying on 
the physics of the problem, we initially derive a corresponding singular model, featured by two 
singular control elements, an arbitrary objective function, spatial wind field as well as equality 
and inequality constraints. Next, we develop a mathematical framework to compute the co-state 
dynamics and the singular controls in matrix form. The mathematical framework is accompanied 
by an algorithm that bypasses the irregularity issue within the co-state system through a recursive 
sequence. Finally, we compare the analytical solutions to the numerical solutions (due to the non- 
reduced model) to show that both coincide perfectly. 

In line with the scope of the present work, we give a review of the previous analytical works 
on aircraft trajectory optimization (those with the vertical plane hypothesis) and explicate the 
major gap in the literature. 

Ardema solved the minimum time aircraft trajectory in the climb phase by the principle of 
singular perturbation, presumably for the first time (see [4]). Particularly, in [4], the minimum 
time problem has been viewed as a singular perturbation which is due to the existence of fast 
and slow variables within the system dynamics. More recently, in [5], an analogous minimum time 
problem has been solved with mass evolution; however, with no relaxed parameters, constraints, or 
wind consideration. The significant outcome of [5] is that the singular minimum time solution can 
be considered a highly accurate solution for the original problem without simplification. Moreover, 
the singular solution (as found in [5]) is in the form of bang-singular-bang (BSB). In [6], the author 
has presented a detailed singular solution for the minimum time problem with wind considerations, 
however, the evolution of mass has been ignored. In [7], the singular minimum fuel aircraft 
trajectory optimization has been solved for the main flight phases, i.e., climb, cruise and descent. 
Specifically, in each phase, the authors have assumed the control as a single input variable together 
with some matching criteria. In [7], the optimal solutions have been sought in the form of BSB for 
all the studied cases. In [8] authors have presented a singular solution for a simplified un-powered 
descent problem with a cost resembling the environmental emission impact. The indirect solutions 
in [8] have been (partially) compared with those from a direct collocation method for the same 
simplified OCP and a very good agreement has been acknowledged. It is worth mentioning that 
the un-powered descent model in [8] includes inequality constraints on the state vector and the 
scalar control (aerodynamic path angle). In [9], the singular maximum range aircraft trajectory 
optimization in the cruise phase has been solved. Specifically, for simplified system dynamics, 
authors have derived the singular extremal by Green’s theorem ( [10] and [11]) which applies to 
a singular OCP with two state variables. In [12] and [13], authors have extended the previous 
work, [9], by accounting for the compressibility effect as well as fixed arrival time in a cruise phase. 
In [14], authors have presented solutions for a 2D aircraft trajectory optimization (vertical plane) 
in the climbing phase with cost being a combination of time and fuel consumption, however, for 
some limited range of the parameters involved and with no wind consideration. More specifically, 
in [14], a direct method has been employed for the initialization of the multipliers associated with 
the PMP model. In this view, the solution procedure in [14] is a hybrid approach. In [15], the 
minimum fuel cruise phase problem with a fixed arrival time in the presence of an average wind 


field is solved. In [15], the optimal structure is sought through BSB for the scalar control, and 
throttle setting. In [16] authors have solved the maximum range un-powered descent problem with 
the aerodynamic path angle being the scalar control under the assumption of BSB optimality. The 
relaxed parameters in [16] are the aircraft’s initial weight as well as the wind parameter. 

In short, all the previous analytical works on aircraft flights (optimization in a vertical plane) 
address the problem with one active singular control on a case-by-case basis and in this respect, 
there is no analytical work considering two active singular controls. In other words, the other 
(possibly) singular control appearing in the problem has been prefixed as a trivially-bang solution 
in the previous works. Here, we show that for a given flight phase, the other singular control 
may switch. Besides, the BSB structure may not always be the optimal structure. Apart from 
the structure of the singular solution, from the available literature, there seems to be no detailed 
comparison between the optimal solutions due to the reduced singular model and those from the 
non-reduced model. This gap is also replenished in this work. 

In principle, the appearance of several singular controls complicates the solution procedure. 
In this respect, analytical works on OCPs with several singular controls are limited in the area of 
optimal control. For this, we address to [17], [18], [19] and [20]. More specifically, in [17] and [18], 
the collision avoidance OCP is tackled analytically. The OCP in [17] and [18] is a singular dual- 
control OCP and authors have shown that the two controls cannot be singular simultaneously. 
Similarly, combination therapy for tumors has been formulated as an OCP with two active singular 
controls in [19] and [20]. In [20], authors have rejected the possibility of a singular surface and 
in [19], authors have used the argument for a more complex model and conjectured, based on 
observations, that the singular surface is not a case of optimality. 

We construct the paper by first defining the original aircraft trajectory optimization model in a 
vertical plane(P“) and the corresponding singular model (P')) in Sec. 2. The novel contributions 
of the present paper have been organized as: 

1) In Sec.3, presenting a compact closed-form solution in form of Jacobian compositions for Po 
with full dynamics, two active singular control elements (aerodynamic path angle and throttle 
setting), spatial wind field, arbitrary objective as well as equality and inequality constraints (the 
closed-form solution can be computed easily in any platform supporting symbolic calculations such 
as MATLAB). 

2) In Sec.3 (Appendix A), detailed insights to the closed-form solution through matrix analysis. 

3) In Sec.4, presenting an effective algorithm to capture various combinations of bang and singular 
extremals and the possible switches inside or outside the singular arcs (the algorithm circumvents 
the irregularity issue within the co-state system). 

4) In Sec.5, solving an example of partial climb phase with full dynamics, equality and inequality 
constraints, wind field as well as two active singular control elements and the objective being a 
combination of arrival time and fuel burn (this cost function is also new in view of the analytical 
aircraft trajectory optimization). 

5) In Subsec.5.2.3, showing that the optimal aerodynamic path angle may follow structures other 
than BSB, i.e., Bang-Singular-Bang-Bang (BSBB) and asymptotic Bang-Bang (aBB). 

6) Subsec.5.2.4, showing that for a pure climb phase, the optimal throttle setting may switch and 
a preset trivial bang structure is not optimal. 

7) In Subsec.5.2.1-5.2.5, comparing the solutions from P'?) with those due to P™). 


2 Problem Statement 


In this section, we introduce the aircraft trajectory optimization problem in a vertical plane 
(original model, P“)). The singular model (reduced model, P')) is presented afterwards. The 
singular model is a reduced model through applying singular perturbation of the zeroth order, [21], 
and order of magnitude analysis in accordance with the physics of a commercial aircraft. 


2.1 Aircraft Trajectory Optimization Problem in a Vertical Plane (Orig- 
inal Model, P“) 


The problem can be formulated as an OCP through the following Bolza form ( [22], [23]): 
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In Eq. (2), x(t) is the tangential coordinate (horizontal distance traveled), h(t) is the normal 
coordinate (altitude), v is the aerodynamic speed, y(t) is the aerodynamic path angle, II(t) is the 
throttle setting, m(t) is the aircraft mass, w is the horizontal wind field, D is the aerodynamic 
drag, T is the thrust force, € is the thrust angle of attack (which is assumed to be zero), L is the 
aerodynamic lift defined as: L = $C; p(h)sv”, C, is the fuel flow, g is the gravitational acceleration, 
C{ is the lift coefficient and t is time. Moreover, Vo4g and Mach number are defined as: 
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The horizontal wind field is a steady-state flow and is a function of coordinates x and h. More- 
over, the cost functional is relaxed, convertible to various forms such as minimum time, maximum 
distance traveled, minimum fuel, as well as minimum environmental impacts or a combination of 
these terms. 

P\) is an OCP with constraints on state and control vectors. In addition, the OCP is singular 
on II(t) and regular on C{(t). 


2.2 Commercial Aircraft Trajectory Optimization Problem in a Vertical 
Plane (Reduced Model, P”?) 


For a commercial aircraft, the significant time scale separation makes it possible to apply the 
singular perturbation technique: 
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In Eq. (5), p is a perturbation parameter. The zeroth term is linked to p = 0. 


Now, we apply order of magnitude analysis to determine the significance of the terms involved. 
For this, we note that Eq. (5) is significant in climb or descent phases. Therefore, it is reasonable 
to take the following orders for the involved variables: sin y = O(107~') (cosy = O(1)), v = O(107), 
Aw = O(10) (w = O(10)), Ax = O(10°), Ah = O(10°) and g = O(1). Moreover, we set € = 0. 

With the above orders, the last term in Eq. (5) is negligible through the following analysis: 
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Therefore, the last term in Eq. (5) can be canceled, leaving: 
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Plugging Eq. (7) into Eq. (2), together with « = 0 and y(t) << 1, P“) turns into the singular 
model, P?): 
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£ is linearly dependent on II if the cost functional contains the minimum fuel term. Therefore, in 
general, this is a singular OCP on both control elements, I(t) and y(t). 


3  Closed-Form Solution for P? 


In this section, we initially present P‘) in vector form. Then, through theorem 1 we show that in 
the case that the state inequality constraints are inactive, a singular surface does not exist. The 
elements of the closed-form solution are specified through a schematic (Fig. 1). In subsection 3.1 
and 3.2 these elements are found in the form of Jacobian compositions. In subsection 3.3 we give 
complementary theorems to deal with the continuity of the co-state variables. 

In vector form, P’) can be recast as: 


min J = L(X,U jdt (10) 
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Where: 
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F = [F\, Fo, F3, Ful” (12) 
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Where, X(t) — [Az(t), An(t), Av(t), Am(t)]", p(t) = MO, MAT, way. war)’, Yo € R* and 
vs € R¢ are the scalar multipliers associated with the boundary conditions. Moreover, 7(t) € R* 
is the multiplier vector for the state inequality constraint. 

The Hamiltonian, associated with Eq. (15), is: 


H(t) = £(X,U) + dP F(X,U) + pC) +n S(X) (16) 


The PMP optimality conditions are: 
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The system is autonomus. Therefore, along the optimality we have: 
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Eq. (19), together with the Hamiltonian condition in Eq. (18), gives: 
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Using Eq. (16) and Eq. (17), the optimality conditions become: 
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In order to obtain y;(t) and II,(t), we re-write the present affine system as: 
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Therefore, the Hamiltonian can be recast as: 


H=L+AMF +. Fy, +A Fan tp C+ S (30) 


Theorem 1. In the case that the state inequality constraints are inactive, t.e., 7 = 0, there is no 
tame interval over which the control elements are both singular. 


Proof. By contradiction, assume that there exists an interval t € |t,,t,] such that the control 
elements are both singular. ‘This means, 
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Plugging Eq. (33) into (32), one ends up with 7 independent algebraic equations for defining 
the optimality with 6 unknowns, i.e., \ € R* plus y,(¢) and II, (t). Therefore, the algebraic system 
is inconsistent and the proof is complete. LI 


Some Definitions: 
1) Singular extremal: an interval over which at least one control element is singular. 
2) Over the singular extremal: U, := P(X, AZ). 
3) Over the singular extremal: AJ := G(X, Az). 
4) Switching point: a time instant at which at least one control element changes its structure. 

There is a possibility that U, has no dependency on A? . In this situation, the solution procedure 
is simplified by removing the loop processing. However, in a more general manner (including the 
present case study in Sec.5), the loop processing is unavoidable. 

The closed-form solution for P) includes finding two significant operators, G and P. With G 
and P found, it is possible to solve over the singular extremal. In this sense, we may suppose that 
the optimal solution contains a singular arc (U,) between two switching points, s) and s(kt+)), 
In this situation, we can solve for t; € [s“*), s+]. Fig. 1 shows this procedure in schematic. In 
the following subsections, we find these operators in closed-forms. 


3.1 y-Singular Scenario 


There is an interval, t © |ta,t)| such that y(t) is singular. Therefore, from theorem 1, H(t) must 
be bang with possible switching instants, i.e., I(t) = Wy := Umaz V Umin. 
For this case, the optimality conditions are: 
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Figure 1: Solution procedure over the singular extremal using the closed-form operators, G and P 


Let us initially assume that the state inequality constrains are inactive, i.e., 7 = 0. Therefore: 
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Plugging Eq. (27) into Eq. (36), y(t) is canceled out and the algebraic system, in Jacobian 
form, defining the co-state dynamics become: 
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The vector, A, is the same as the commutativity of the vector functions in Lie bracket notation 
(see [14]): 
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Eq. (45) can be expanded as: 
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The time derivatives of of. and A are: 
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Eq. (46), together with Eq. (27) and Eq. (47), gives: 
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For a non-saturating singular control, the generalized Clebsch (Legendre) condition ( [24]) for 
the present OCP is equivalent to: 
0 d?S7 
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In the case that the state inequality constraints are active, i.e., 7 # 0, it can be assumed an 
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Since, we are necessarily in the singular arc territory, the following holds: 
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For t € |tc,ta] we have, ys = Ys,a-. Therefore, 7, can be determined in terms of the state 
variables and A, through the time derivative of Eq. (54). However, this derivation is omitted 
since the value of 7 is out of necessity and 7¥s,_ suffices to process the integration. 
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3.2 Il-Singular Scenario 


There is an interval, t € [ta,ty| such that I(t) is singular. Therefore, from theorem 1, y(t) must 
be bang with possible switching instants, i.e., y(t) = Yo := max V Ymin- 
For this case, the optimality conditions are: 
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We initially assume that the state inequality constrains are inactive, i.e., 7 = 0. Therefore, 
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The solution of the co-state system, Eq. (57-59), for An, Ay and A», reads: 
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Eq. (65), together with Eq. (56) and Eq. (66), gives: 
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The generalized Clebsch (Legendre) condition requires: 


ie) 


_ 288" __yr(yp)F. 1 — JF. q)B) —TL>0 (70) 
dll dt2 — sll sll ee 
In the case that the state inequality constraints are active, i.e., 7 # 0, it can be assumed an 
interval, |t/.,t)| C [t,,t,] such that at least on element of S is zero; i.e., Sg = 0. Therefore, for 


t € [t, &|, from a = il), 


(St) (F + ws.) 
(Sx) Fon 


Nq is connected to the co-state dynamics through: 


II,a(t) = — we a (71) 


OLsun/z ” 
_ MB + a 4 + pF ayy) — a Fam 
Na = dS. 
OX s,Il 


(72) 


3.3. Jump Conditions and the Complementary Theorems 


The co-state variables may be discontinous at any time 7 in the boundary arc interval [25]. There- 
fore, the jump conditions for any contact time 7 read, 


OS 


Mie Ve 2 ae ha u(r) ae v(t) >0, w'(r)S=0 (73) 


Since S = S(h,v), for the active element of S, i.e., Sg = 0, Eq. (73) gives, 


OSa 


ah (74) 


An(t~) = An(7*) — u(r) 


OSa 


Av(T~) = Av(TT) — v(7) Fn (75) 
Ae SAT) (76) 
Mae Sena”) (77) 


Theorem 2. (Junction theorem in [26]) Let ta be a time that the singular and nonsingular sub- 
arcs of an optimal control u are joined. Let q be the order of the singular arc. Suppose the 
strengthened GLC condition is satisfied at tg, and assume that the control is piece-wise analytic 
in the neighborhood of t,. Let u\") be the lowest-order derivative of u that is discontinuous at tg. 
Then g+r is an odd integer. 


For the present OCP, q = 1. In addition, due to the bang/singular nature of the optimality, 
we suppose that the control vector itself is discontinuous, i.e., r = 0 for singular arc + nonsingular 
arc interaction. Therefore, g +r = 1, is an odd integer. 


Theorem 3. (theorem 5.1 in [27]) Let ty be a time that an interior non-singular arc and a 
boundary arc of an optimal control u are joined. Let u\") be the lowest-order derivative of u, that 
1s discontinuous at ty, and let p be the order of the state inequality constraints. By the assumption 
ofp <2q+r, if v(ty) > 0, then p+r is an even integer. 


Theorem 4. (theorem 5.4 in [27]) Let t. be a time that an interior singular arc and a boundary arc 
of an optimal control u are joined. Suppose that the strengthened GLC condition holds. Let u“, 
be the lowest order derivative of u, that is discontinuous at t.. By the assumption of p < 2q+r7r, 
V(t.) =0 andq+r is an odd integer. 


For the present OCP, p = 1; hence, p < 2q +7 is always satisfied. The boundary arc remains 
in the interior of the admissible control set and in a bang/singular optimality, it can be figured as 
a part of the singular arc switching between |y,,II,]* and [ys5.a, s.a]*. Therefore, it is reasonable 
to assume that the control vector itself is discontinuous in the nonsingular arc+-boundary arc and 
singular arc+boundary arc interactions; so, r = 0. From theorem 3, since p+ r= 1 is an odd 
integer, so v(t,) = 0. Theorem 4 is automatically satisfied and so, v(t.) = 0. Therefore, we can 
conclude that the co-state variables are continuous along the optimal trajectory. 


3.4 Simplified Cases of the Closed-Form Solution 


For some special cases, it is possible to extract the singular arc expressions from the closed-form 
solution. Through theorem 5, 6 and 7, we show that the singular arc formulas are linked to 
det|M] = 0 or det|M’] = 0 (see Appendix A). 

Extraction of the singular arc formula is of help to predict the instants of the junction points 
in a presumed optimal structure. However, as shown in Appendix A, it is doable only for some 
special cases. Therefore, in principle, the presented closed-form solution still lacks identification of 
the number of the junction points and the corresponding instants, pending for a closing theory in 
optimal control. Loosely speaking, a shooting procedure is typical of this situation. In general, the 
singular control can be assumed as, U, = P(X, A*). However, for some cases, this form is reduced 
to U, = P(X). The reduced form, i.e., U; = P(X), leads to a considerable simplification in 
devising a solution algorithm. It is worth mentioning that the previous cited works on analytical 
aircraft trajectory optimization in a vertical plane fall within the category of the singular arc 
formula or U, = P(X). 

In the following section (Sec.4), we consider the case of U, = P(X, *) which is new and more 
general, and also it includes the case study that we have solved in the present work. 


4 Solution Approach for P’) 


A multiple shooting procedure for the present OCP can be figured as: 
O(Y) = (sur any Onail’, OE Rt 


a - . 78 
Y S(t, ate rAe@)|, inate YeR’™ (78) 


13 


In Eq. (78), Ox,k = 1,...,N +1 are the boundary (closing) equations and t,,k = 1,...,.N are 
the number of junction points. In addition, t* is a junction point joining a bang extremal to a 
singular extremal. 

In a conventional manner, the closing equations are, 


O1.q<a = bf (Xp) = Xp —X(tp) = 0 
Og41 = H(ty) = 0 


In the case that q # 4, we can consider A4(t) = ... = Agsil(tz) = 0. Therefore, normally, 
© € R°. This means that a standard shooting procedure can handle a limited number of the 
shooting parameters, i.e, Y € R°. However, there may exist cases with additional shooting 
parameters, e.g., an optimal solution with BSB structure for both y(t) and II(t). 

Therefore, in general, in order to automate the solution procedure (to extricate it from searching 
for amenable closing equations), we propose an alternative solution approach; that is to translate 
the shooting problem into a very sparse parameter optimization in the following form: 


_min TI = &(X (tx), tx) 
Aa (t*), ty, k=1,..., KEN 


(79) 


S.t. 


dX . 80 
= F(X, Ao(t*)) oy 


bo(X (to), to) = 0 
7 (X(tk), tx) =0 
First, we interpret the representation of the junction points in a parameter optimization. Next, 
for a minimum time problem, we compare the shooting results with those from a corresponding 
parameter optimization and draw the conclusion that, conversion of the shooting problem to a 
very sparse parameter optimization surpasses in terms of labor, automation and also insertion of 


the boundary conditions for y(t) with an analogous level of accuracy and computational time (see 


Appendix A). 


4.1 The Proposed Algorithm for P) 


With the help of theorem 1, we design an algorithm (see Algorithm 1) that the unknown parameters 
are only the junction points, t, and \,(t*). Algorithm 1 allows the possible switches inside and 
outside the singular extremals. The maximum number of singular extremals are assumed to be 
one for y(t) and one for II(t). Algorithm 1 (in its current form) is applicable to the simulation 
of a mixed climb/cruise phase. However, with slight modifications, the logic is extendable to the 
case of mixed cruise/descent phase and free routing flights. 

With A* known over the entire time domain, it is straightforward to verify the optimality of 
the singular extremals by computing Eq. (51) and Eq. (70). In addition, the switching functions 
can be computed to check the first order optimality condition of the bang extremals. Initializing 
the bang extremals for yy2:)5 and IIy1.,4 seems to be a try and error procedure and in most cases, 
we verified that yo =... = ypx and II,; =... = I[y4. However, in practice, we were able to obtain 
the optimal solutions by initiating the algorithm with different combinations of the extremals. In 
this situation, we sometimes witnessed some pulses in the optimal solutions which were linked 
to the wrong extremals. More specifically, suppose that the preset 7; over (tj41 : ti:2) is not 
optimal. In this case, the optimization module computes (t;41 & tj2) secured to some decimal 
places. This very marginal deviation shows itself as a pulse. Nonetheless, a verified structure 
preserves itself in a range of variations over the relaxed parameters involved. 


5 Case Study 


In this section, we set the objective functional as to be a combination of the arrival time and fuel 
burn and solve the OCP for a partial climb phase. In the following, we explore the appearance of 
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Algorithm 1 The Multi-Structure Detection Algorithm 
Figuring the extremals as: 
(0,1) : Ye1, oi 


(t1,t2) : Ys, Ibe 
(t2,t3) : Ys, [lbs 
(t3,t4) : Yeo, Hoa 
ee 73, ITs 


cK 

On 
ch 

on) 


Il, 

(tg, t7 := ty) : Yos, [os - - 

The decision variables are: th :k = 1: 7,Az2(t1). 

Assumption: 

Ys,a = Ys 

Hea S IT, 

Initializing the bang extremals: 

Veist =1:9¢ Vein V i) 

ae) —-1:5¢ (Donia V Miiae) 

Initializing the decision variables: 

fy = BOR =1:7,r0(E1) = AL? (EO) 

Integration process: 

a) Solving 4* = F(X,U) over (0, f°) 

b) Solving 4* = F(X,U) and one = aN simultaneously over (£0), #) with: 
Nis Noss Ney, = G(X Xe) 

ye = PO(X,2) 

c) Solving a = F(X,U) by forward integration from £9) to i) and dx = a by backward 


» 
oa 
iS 


M 
M 


and forward integration from f°) to £0) =: 0 and #9) to i with: 
1, 2= P(X) 
Solving the associated NLP 


the various scenarios captured by Algorithm 1. More specifically, by changing a from 1 to 0, the 
OCP alters from the pure minimum time to the pure minimum fuel. For the simulated range of 
the parameters involved, a = 1 corresponds to a uniform B structure for II(t), i-e., II*(¢) = Uma. 
However, the optimal structure for y(t) includes one singular extremal. It is noteworthy mentioning 
that a = 1 corresponds to the case of one active singular control (as mentioned before, the case of 
one active singular control has been the target of numerous previous works). As a changes from 
1 to 0, the optimal I(t) can be BB or BSB; whilst, the optimal y(t) is mostly BSB. 

For the comparison purpose, the optimal solutions due to PMP are labeled as the analytical 
solutions and those from the direct approach are labeled as the numerical solutions. 


5.1 Problem Setup 


In order to engage the full dynamics, we simulate a partial climb phase. ‘The objective is considered 
as a combination of the flight time and the fuel burn: 


L=a+(1-a)NIC.T, a<0<1 (81) 


T, C, and D are simulated by Base of Aircraft DAta (BADA) model from EUROCONTROL 
providing a general smooth aircraft performance model [28]. In addition, air density is approxi- 
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mated by International Standard Atmospheric (ISA) model: 


T(h) = Cr, (1- = +h?Cr,), Cs(v) =Cs,(1+ a ) 
2 P(S)™, O(h) = Oo — Bh, p(h) = mah (82) 


1 
Di, h). = 5 oth) su" (Cp, + Cp,C7) 


Where, s is the aerodynamic lift surface and p is the air density. C'p,,2 = 1, 2,3, s, Cp,,7 = 1, 2, 
C,,, = 1,2, R, 8, Po and 0 are known constants. 

The model constants as well as the fixed boundaries for different scenarios are tabulated for a 
medium-haul aircraft (see Table 1). 


Table 1: The model constants and the fixed boundary conditions 


Parameter Value Unit (SI) Parameter Value Unit (SI) 
h(0) 3480 m v(0) 151.67 m/s 
h(t) 9144 m v(t pf) 191 m/s 

8 122.6 m? g 9.81 m/s? 
Cr, 141040 N Cz 14909.9 m 
Cr, 6.997e~ 1° m* Cp, 0.0242 — 
Cp, 0.0469 — Cs, 1.055e-°  kg.s-+ N71 
Cas 441.54 m.s* R 287.058 J.kg7!K7! 
Qo 288.15 K 8 0.0065 Km! 
Po 101325 Pa x(0) 0 m 
Mer 0.8 — mit) relaxed kg 


We consider a rather arbitrary wind field as: 


WG, hy—G,,(0,,h° G,, X10 “2° 410-2) (83) 

The adjustable parameters a®,,al, and a?, can produce different wind profiles ranging from 

smooth to severe pure positive (or pure negative) winds as well as mixed positive/negative winds 

of smooth or severe types. The latter case is to naively imitate the presence of storms. Some of 
the distinct profiles simulated in this work are graphed (see Fig. 2). 


w(a,h):a°, = 1, a), = 0.2 x 107°, a2, = 0.07 w(a,y) : a®, = 2,a!,=0.2 x 10-6, a? = 0.2 
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Figure 2: a view on two different wind profiles duo to Eq. (83) 
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5.2 Results 


The results of the present simulation are given through subsec.5.2.2-5.2.5. The direct numerical 
solution procedure for P“”) is specified in subsec.5.2.1. 


5.2.1 Numerical Solution for P 


For the numerical solution, we solve the complete model (sec.2) based on Adams-Bashforth [29] 
discrete scheme of the second order with a piece-wise constant assumption for the control vector, 
U = [C\(t), 1 (t)|’. The resulting NLP can be figured as: 


= at b=) Xie) Sa 
ee On <6 ( a) ( a 0) 4( Nn)) 
S.t. 


dX 
a F(X,C,,ll), X,F ER? 


®o(X(to)) =0, GER? 
®v(X(tv)) =0, On € R’ 

Of X5 <oX5 mae, Unin SU <1 
G(X)<0, GER 


(84) 


Where, X = [z,h,v,m,7|’. In order to simulate the partial climb phase, the boundary 
conditions for Xs are considered as: 


X5(to) = X5(tnw) = 0 (85) 


For the numerical solution, interior-point /barrier algorithm of the NLP solver fmincon from 
MATLAB@®) Optimization Toolbox, was adopted as the optimization module. 


y(t) H(¢) 


——N ON 
5.2.2 Scenario 1: BSB+ B 


In this subsection, we explore BS B+ B structure that is the case with one active singular control. 
We consider the case with a = 1 corresponding to the minimum time problem. For the simulated 
range of the parameters involved, the optimal solution for II(t) was confirmed to be a uniform bang 
structure; ie., I*(¢) = I, = mae = 1. For the case of NoWind (a?, = 0) with m(0) = 69000, 
Ymin = 0, Ymax = 0.2 and varying x(t), the analytical solutions for y(t) and the state variables 
are compared with the numerical ones in Fig. 3. From these plots, it becomes clear that the 
singular and the boundary extremals are perfectly covered. 

Fig. (4) shows a comparison for y(t) between the analytical solutions and the numerical ones 
for two distinct cases with wind consideration. 

A better comparison is to compare the optimal objective and the final mass. In this respect, 


these differences between the analytical and numerical solutions are tabulated (see table 2). The 


= aaa = pene Am | = ail = ia meee 


differences are defined as, |At | 
Interestingly, the deviation is negligible. 


y(t) II(t) 


———_——_s 
5.2.3 Scenario 2: BSBB,aBB+ B 


In this scenario, a = 1 and for the simulated range of parameters, II*(t) was confirmed to be a 
uniform bang structure; i.e., II*(¢) = I, = 1. The BSBB structure may exist for a range of the 
involved parameters. However, in the present simulation, BSBB for 7(t) was captured in larger 
values of Ymin and specified mr. AS Ymin increases, the optimal structure for y(t) changes from 
BSB to BSBB. Thereafter, the new structure, BSBB, is dominant and the singular extremal shrinks 
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Figure 3: 7(t) and the state variables: the impact of x(ty), a = 1, NoWind (a¥, = 0); Red Lines: 
Analytical Solutions, Black Lines: Numerical Solutions 


Table 2: |At;| and |Am,| for BSB+B structure; mp = 69000, Ymin = 0, ¥max = 0.2,a9, = 0 
zy |Atyl(s)__ |Amayl(&g) 
140000 122 0.27 
160000 1.17 0.28 
180000 Lo 0.34 
200000 1.42 0.32 
220000 1.26 0.35 
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mo = 69000, ymin = 0; Imax = 0.2, 21f = 180000 


0 100 200 300 400 500 600 700 
t(seconds) 


(a) 


Figure 4: y(t): the impact of wind, a = 1, (ai, = 0.2 x 107°, a?, = 0.07); Red Lines: Analytical 
Solutions, Black Lines: Numerical Solutions 


aS Ymin increases. ‘This shrinkage continues until the singular extremal disappears asymptotically; 
turning BSBB to aBB in a critical ymin, 1-€., Ymin,er.. Therefore, aBB is the closing optimal 
structure since for min > Ymin,cr., the number of boundary equations are more than the shooting 
parameters and therefore, in principle, the shooting system is inconsistent. 

For the sake of brevity, we only graph the main variable, y(t). Fig. 5 shows a comparison 
between the analytical solutions and the numerical ones for fixed ym; and a. = 0 (Fig. 5-a) 
and varying ymin (Fig. 5-b). Fig. 6-a and Fig. 6-b display the analytical solutions for positive 
and negative wind profiles respectively with Ymarz = 0.15. Fig. 6-c shows the analytical solutions 
indicating the evolution of BSBB to aBB for a?, = 0. In addition, table 3 depicts |At;| for some 
stages with mo = 69000, Ymaz = 0.15, mr = 68100, x¢ = 140000. 


mo = 69000, x ¢ = 140000, ymin = 0.03, WindIncluded 


(1) : a®, = 0.5; (2) : a®, = -0.5 


mo = 69000, x = 140000, ymin = 0.03, 0.032, NoWind 


0 100 200 300 400 500 600 0 100 200 300 400 500 600 700 
t(seconds) t(seconds) 


(a) (b) 


Figure 5: y(t): BSBB comparison; Red Lines: Analytical Solutions, Black Lines: Numerical 
Solutions (my = 68100, at, = 0.2 x 10~°, a?, = 0.07) 


Table 3: |At;| BSBB+B structure; a?, = —0.5, aj, = 0.2 x 10-°, a2, = 0.07 
Ymin |At¢|(s) 


0.027 0.93 
0.029 0.95 
0.031 0.91 
0.033 0.97 
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0.14 mo = 69000, x; = 140000, a®, = 0.5, al, = 0.2 x 10°°, a2, = 0.07 0.14 my = 69000, z+ = 140000, a;, = —0.5, ay, = 0.2 x 107, ai, = 0.07 


Ymin = 0.025 : 0.01 : 0.032, Ymazx = 0.15 z 
— 


Ymin = 0.025 : 0.01 : 0.032, Ymax 
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mo = 69000, m¢ = 68100, ymin = 0.027 : 0.01 : 0.036, 0.0365, 0.03668, NoWind 
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Figure 6: y(t): BSBB; m+ = 68100, aj, = 0.2 x 10~°, a2, = 0.07 
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5.2.4 Scenario 3: BSB+BB 


In order for BB structure of II(t) to appear, the involved parameters need to change. In this 
respect, we set a = 0.7 and adjust the other involved parameters to emphasize the existence of BB 
structure for II(t). For the case that the state inequality constraints are inactive, i.e., 7(t) = 0, 
the optimal controls as well as the corresponding optimal states are graphed in Fig. 7. ‘The 
effect of ILnin with active state inequality constraints and qa = 0 is shown in Fig. 8. Fig. 9 
shows the impact of wind on the evolution of the optimal solutions with II,,;, = 0.9. Fig. 10 
displays a comparison between the analytical solutions and the numerical solutions for the case of 
ILnin = 0.9. Table 5 shows the results for |Aty| and |Ams| clarifying a remarkable resemblance 
between the two solution approaches. 


Table 4: |At,-| and |Am,| for BSB+BB structure; mo = 59000, xf = 150000, min = 0,%max = 
015.0), 1.4, = 02% 10-°,4 = 0.07 
Tin [Atsl(s)__ |Amyl(kg) 


0.95 1.12 0.15 
0.9 0.98 0.13 
0.85 0.93 0.11 
0.8 1.03 0.13 
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mo = 59000, x¢ = 150000, aj, = 0.2 x 10-°, a?, = 0.07 
a= Inn = 0.7 
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Figure 7: BSB+BB; n(t) = 0, mo = 59000, Ymin = 0,¥max = 0-15, Umin = 0.7, 2¢ = 150000, a}, = 
0.2% 10° a7 = 0.07 


y(t) II(t) 
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5.2.5 Scenario 4: BSB+BSB 


In this scenario, we set a < 0.25 and, in particular, examine the effects of I,in, a and wind. Fig. 
11 shows the impacts of II,,;,, and a@ in the presence of wind. From these plots, as a decreases, the 
singular extremal for I(t) expands. Fig. 12 shows the impact of a severe mixed positive/negative 
wind with IL,jn = 0.3 on the evolution of BSB+BSB structure. The wind profile used in the 
production of Fig. 12 naively captures some of the characteristics of a storm. This, by itself, 
highlights the applicability of the proposed analytical approach to handle complex phenomena. 
In this scenario, the numerical solutions show severe chattering response due to the fact that 
the OCP is singular on I(t). We only show a comparison between the analytical solutions and 
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mo = 65000, x + = 170000, NoWind, a = 0.7, Me = 0.8 


y(t) (rad.) 


I 
0 100 200 300 400 500 600 700 “0 100 200 300 400 500 600 700 
t(seconds) t(seconds) 


(a) (b) 


mo = 65000, x» = 170000, NoWind, a = 0.7, M., = 0.8 


y(t) (rad.) 


0.05 


0 100 200 300 400 500 600 700 


(c) (d) 


Figure 8: BSB+BB; mo = 65000, ymin = 0, ¥max = 0-15, x77 = 170000, a2, = 0: (a) and (b) belong 
to min : 0.95 — 0.05 : 0.7. (c) and (d) belong to ILnin : 0.6 — 0.1 : 0 


the numerical ones for y(t) and II(t) and this suffices to clarify the superiority of the analytical 
solutions (see Fig. 13). Table 5.2.5 is devoted to a comparison for |Atr| and |Am,|. Unexpectedly, 
the deviation is clearer. This deviation is perhaps due to the inaccurate results from the numerical 
solutions. 


Table 5: |Atr| and |Am,| for BSB+BSB structure; mo = 59000, 77 = 150000, min = 0, %max = 
(0.15 @).= le, =02 <10-",a, = 0.07 
a,Umin  |Aty|(s) — |Ams|(kg) 


0 7.68 0.86 
0.1 4.32 0.67 
0.25 2.93 0.45 


6 Conclusion 


It was presented a mathematical framework for commercial aircraft trajectory optimization in 
a vertical plane with two active singular controls, an arbitrary objective function, equality and 
inequality constraints as well as arbitrary wind profiles (Sec. 3). Conventionally, a shooting proce- 
dure is typical of this problem (Po. However, we translated the shooting procedure into a very 
sparse NLP and gave mathematical justifications for the correspondence (Sec.4). An algorithm 
(Algorithm 1) was proposed, able to capture multiple combinations of bang/singular extremals. 
Algorithm 1 is also able to capture the effects of boundary conditions for y(t) through a zeroth 
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Figure 9: BSB+BB; mo = 59000, ymin = 0, maz = 0.15, x ¢ = 150000 


order expansion scheme which includes additional junction points. Specifically, we compared the 
optimal results due to (P)) with those from (P“)) with the overall conclusion of remarkable 
accuracy. Through the case study, we showed in Sec. 5.2.4 that a pure climb phase (not ending 
up with a cruise phase) may involve optimal II(t) in the form of BB with switching instants inside 
or outside the singular y(t). More specifically, this new structure can be understood through the 
objective function which was a combination of the arrival time and the fuel burn. In addition, we 
showed in Sec. 5.2.3 the existence of optimal structures for y(t) other that BSB, i.e., BSBB and 
aBB which were also new in view of the analytical aircraft trajectory optimization. 

It is well-known that the optimal solutions by PMP involve some complications especially in 
complex problems and inherently, PMP approach to OCP is more difficult, in terms of labor 
and implementation, than a direct NLP approach. A good deal of this difficulty is addressed 
to the initialization of the multipliers (whether to be the state equality multipliers, co-states, or 
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Figure 10: BSB+BB; mo = 59000, ymin = 0,7 max = 0.15, 75 = 150000 


those from the inequality constraints). Nevertheless, for a singular OCP, the PMP approach may 
lead to a considerable simplification. On the other hand, in a direct optimization, the projected 
Hessian matrices are ill-posed if the problem is singular. Hence, concerning the accuracy and the 
computational time for a singular OCP, the PMP approach prevails. 

The presented closed-form solution for P‘) was verified to be quite efficient and accurate 
especially when the optimal structure is identified. The notable point is that an identified structure 
likely preserves itself over a range of the relaxed parameters. It is hopeful that the present work is 
evident enough that the singular PMP approach is applicable to more complex arenas regarding 
the commercial aircraft trajectory optimization. 
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Figure 11: BSB+BSB: the effects of [min and a: mo = 59000, ymin = 9, max = 0-15, x = 150000 


A Appendix A 


Theorem 5. /f \,(t) = 0, the singular arc of y(t), is defined by det|M] = 0 and the occurrence 
is subject to L=c(az,m) Fy (F:= F+0F, 1). 


Proof. For the singular arc to be defined purely in terms of the state variables, we rewrite the 
algebraic equations (Eq. (37),(38), (39)), taking into account, A,(t) = 0 and Fyiy = Fay = Fo = 
OQ, as: 


Ah Ag F's3, A3Fy — AgF 3 —FiF 33,4 —£ 
det |M] Av — Ail oa —A2F4 Py F324 a (86) 
Am A3 Fis2,4 — AoF’s3, Aoks — £3 F597 (ae) Fs, 
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Figure 13: BSB+BSB: comparison between the analytical solutions (red lines) and numerical 
solutions (black lines) 


For det|M] to be zero, the right hand side of Eq. (86) needs to vanish. Therefore: 


<~) Fey =0 (87) 


0 (88) 


ee =) (89) 
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On the other hand: 
det|M/] = A3FyFs9,y = Ag Fy F.3, = Aa Fs Foo, (90) 


The vector elements of A are: 


“w 


jae OF; 
A; — eae ono fh _ pa ae ay Pein = i, 2, 3,4 (91) 


Eq. (87), or equivalently Eq. (88), together with A,, from Eq. (91), give: 


“w~ 


L=cean)k (92) 


With Eq. (92), it is straightforward to verify that Eq. (89) corresponds to det|M] = 0 and the 
proof is complete. L 


Corollary 1. For the minimum fuel OCP, L:= —Fy. Therefore, from theorem 5, the singular arc 
for y(t) (I(t) = Il, ) can be defined in terms of the state variables by det|M] = 0 if Ax (t) = 0. This 
condition corresponds to the case of unconstrained horizon, i.e., unspecified x(tr) and w = w(h). 


Theorem 6. For the objective function L:= —F, = —IIF 4, the singular arc of I(t) is defined 
in terms of the state variables by det|M"| = 0 if Ax (t) = 0. 


Proof. On using the expressions for M’ and R’ (given by Eq. (62) and Eq. (63) respectively) 
together with Fy = 0 (F := F + %F.,y), Lea = 0 and Az (t) = 0 we have: 


Xp, BaF 33.0 _ b3F sam —ByP F3Fea.n =n 
det/M") Ay = BoF sa BaF» —FoF ean oF aa A (93) 
a a a s,II 
Am —BoFs3 0 Bo F3 — B3F5 FoF 53 11 ax 's,0 — Ox fr 


On setting £ = 0, Leu = —Fsan, the right hand side of Eq. (93) becomes: 


OF 54 WA 

47 aX (94) 
OF 4 8 ee 
aX T+ PA (95) 
x n x OF sa. 5 
Poa tt (BoF3 — B3F2) + FoF 53,0 FX F=0 (96) 
On the other hand: 

det(M"] = BoF3Foan — B3F2F san + BaFoFs3.n (97) 

OF sia 2 OF, 
SS es }= 1,2,3,4 (98) 

B; t=1 OX; 7=1 OX; Il, J ) ,3, 

On accounting 6, from Eq. (98), it is straightforward to approve that Eq. (94) and Eq. (95) 
are automatically satisfied and Eq. (96) is equivalent to det|M’] = 0. O 


Theorem 7. In the case that F, = 0, the singular arc for y(t) is defined in terms of the state 
variables through det|M| = 0 if £ = d(a)F\. 


va 


Proof. Since Fy, = 0, am. = 0. Therefore, the co-state dynamics, Eq. (37), Eq. (38) and Eq. (39) 
become: 


Ne FLO PF; ey 
det|[M] [An] =adj| 0 Fy Fs3,7 0 (99) 
Np Aj Ao» A3 of Fs 


In order to have det|M/|] = 0, RHS of Eq. (99) must vanish: 


n OL 
L( Ao Fs, = A3F's2,7) —_ F3 F595 ax sy = 0 (100) 
~ OL 
= — fF, —F,, = 101 
AiL lays 0 (101) 
~ OL 
a 102 
AiL cia Soil Ox sy 0 ( ) 
On the other hand: 
det|M] = A3FFs2, — Ag FF 53,7 — A, F3F yo, (103) 
and: 
OF 51 A OF, OF, 
=> 2 oe eg 104 
A = 1 OX. t=! OX; Y Ox y ( 


From Eq. (104), together with Eq. (101), or equivalently Eq. (102), 


“w 


L=der (105) 
Plugging Eq. (105) into Eq. (100), gives det|M] = 0. O 
Corollary 2. On setting d(x) = —1 in theorem 7, we have, £L = —F, = —(u+w). On the other 


hand, Fy, = 0 is equivalent to I(t) = 0. Therefore, it becomes lucid that this case is applicable to 
an un-powered descent with the maximum range, being the objective function. 


A Appendix B 


In view of a parameter optimization we recast the shooting problem as: 


_min 7 = &(X (tx), tx) 
Ax (t*), ty, k=1,..., KEN 
S.t. 
AX 
ae F(X, Az(t")) (106) 


do ( X (to), to) = 0 
os(X (tx), tx) =0 


Since we have the expressions for the optimal controls over the extremals, the parameter 
optimization contains only the missing parameters, i.e., A,(t*), and the number of the extremals, 
th. 

The time domain is accumulation of the time intervals: 


K 


0 = |) 2% (te-1, tx) (107) 
k=1 
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The numerical integration is performed over each subdomain, 0, (with an equidistant mesh): 


~ ~ k-1 k 
eee 
Ai, = = =tj-tj1, J=1+5 9 Nj,..,5 NM, b=1,2,.,K (108) 
al ak 


For a better understanding and also to simulate the presence of the state multipliers, co-states, 
we form the Lagrangian as: 


N K 
J = ®(X(tk), tk) +p dotvRek t+ > AMA, N=S My (109) 
j=] k= 1 


A is the state deficiency in the discrete points, t;, verifying the presence of co-states, A: 


~ 


A, = X(tj) — X(t) (110) 


X(t;) is the expanded field about X(t;—1). In order to facilitate understanding the optimization 
process, we take a basic Euler expansion: 


~ 


Kak Gara) ONG =b) (111) 
Therefore: 


J = O(X (tk), tk) + bo + VK OK+ 
ear Ni 


~ T oh ae (112) 
So SO NP ( X(ty-1) — X(t) + F (X (t5-1), Av (#*)) Ate 
R=] fait yy Mi 
The KKT condition for t, instates: 
OF 
— = Q, a eaten 2 © 113 
ai, (113) 


For the interior points, i.e., k = 1,.., A —1, the KKT condition involves the adjacent extremals, 
At, and At;,.1. Therefore: 


eee De (te Aa) TUN Gag et )) 0 
Oty Nz ys J ( (¢; 1); ( )) Nut » 5 ( (t, i) ( )) 
j=14+ D2 Ni j=14+>D"_, Ni 


(114) 


On the hand, the discrete Hamiltonian for this system is H; = A; F(X (t;), Ax(t*)). Therefore, 
on defining a deficient Hamiltonian, H; = Ni; F(X (tj-1), Ac (t*)), one arrives at: 


Ay = Ansys (115) 
With: 
= 1 sai Ni 
H, = a APC), Sate at (116) 
a 
For the boundary point, k = Kk: 
OF O08  ,O0bK , 1 ma T : OD  7(O0bK , Zz 
a a ea SAP F(X(tj-1), Ac(t*)) eo +Hx =0 


Flt NN; 
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In short, Eq. 


(115) and Eq. (116) correspond to the uniformity of Hamiltonian and the 


transversality condition for an autonomous system in Mayer form. 


In the case that N, = 1,k = 1,2,...,.kK = N, ... gives: 


Fea Hag, F212 ayN (118) 


Together with the assumption, N — co: 


~ 


d(H) =0 (119) 


Therefore, it becomes perspicuous that the junction points correspond to the uniformity of 


Hamiltonian over the adjacent extremals. 


It turns out that a very sparse parameter optimization over the parameters, ty, Az (t*) is re- 


markably accurate compared to a shooting translation and hence, it is a logical alternative (see 
table 6). The data in table 6 are from a simulation with £ = 1 and No Wind condition. The 
optimal y(t) and I(t) were confirmed BSB and B respectively. As for the shooting and parameter 
optimization, we used fsolve and fmincon from MATLAB respectively. 


Table 6: Shooting (s.) vs. 


Parameter Optimization (p.o.) for a Minimum Time Problem: No 


Wind, mo=79000, x,¢ = 180000, BSB (parameters: ¢1, ta, t3 := tr, Ax(t1)) 


Quantity Test Case 1 ‘Test Case 2 ‘Test Case 3 Test Case 4 
N, 100 300 500 1000 
No 500 500 500 1000 
N3 500 500 500 1000 
t1(s.) 79.1012 78.1347 77.9744 77.7547 
t1(p.o.) 79.0602 78.1651 77.9946 77.7574 
to(s.) 819.0109 821.0051 821.4188 821.8351 
to(p.o.) 818.9088 820.9962 821.4257 821.8334 
tr(s.) 840.5155 842.2944 842.6612 843.0167 
t (p.o.) 840.5131 842.2944 842.6612 843.0167 
Ar(ti)(s.)  -0.001533 = -0.001246 ~—- 0.001191 —--0.001168 
Nx (t1)(p.o.) -0.001656 = -0.001268 + ~—-0.001180 —--0.001181 
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